#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <omp.h>
#include "common.h"

// This will be the array containing all of the boxes/simulation domains
box ***boxes;
// The root box
box world;

#define THREAD 16
//
//  benchmarking program
//
int main( int argc, char **argv )
{    
	if( find_option( argc, argv, "-h" ) >= 0 )
	{
		printf( "Options:\n" );
		printf( "-h to see this help\n" );
		printf( "-n <int> to set number of particles\n" );
		printf( "-o <filename> to specify the output file name\n" );
		return 0;
	}

	int n = read_int( argc, argv, "-n", 1000 );
	char *savename = read_string( argc, argv, "-o", NULL );

	FILE *fsave = savename ? fopen( savename, "w" ) : NULL;

	particle_t *particles = (particle_t*) malloc( n * sizeof(particle_t) );
	set_size( n );
	init_particles( n, particles );

	//-----------------------------------------------------

	// Creates the main box, this is box number 0, the root
	init_world(&world);

	// Adds the particles to the world box
	for(int i = 0; i < n; i++ )
		world.add_new(i);

	// Calculates what the depth level should be
	// The depth level is the number of layers our simulation will have
	int depth_level = ceil(log10(n));
	printf("depth_level = %d \n", depth_level);

	// Each layer is a 2 dimensional array containing boxes
	// layer 0 has 1 box
	// layer 1 has 4 boxes
	// layer 2 has 8 boxes
	// etc...
	boxes = new box**[depth_level + 1];	// Creates the level 0 multi dimensional array
	
	// Initializes the 2 dimensional arrays, creates the boxes
	initializations(boxes, depth_level);
	// Links the root boxes 4 children to it and sends them their particles
	split(world, particles, boxes, 0, 0);
	// Links the 4 level 1 children with their own children
	// and finishes linking all the other boxes with their children
	// also, finishes distributing all the particles to each box accordingly to their coordinates
	build_tree(boxes, particles, 1, depth_level, 0, 0, 1, 1);

	//
	//  simulate a number of time steps
	//
	double simulation_time = read_timer( );

	#pragma omp parallel num_threads(THREAD)
	{
		// This part will calculate which of the boxes this thread is responsible for
		// Each thread only calculates for a part of the simulation domain
		// This is done at level 2 where there are 16 boxes
		// at maximum this means we can have 16 threads
		int thread_id = omp_get_thread_num();
		int x = thread_id / 2;
		int y = thread_id % 2;

		int avg = 16 / THREAD;
		int mod = 16 % THREAD;

		int index;
		int next_index;

		if(mod == 0)
		{
			index = thread_id * avg;
			next_index = (thread_id + 1) * avg;
		}else if(thread_id < mod)
		{
			index = thread_id * (avg + 1);
			next_index = (thread_id + 1) * (avg + 1);
		} else
		{
			index = mod * (avg + 1) + (thread_id - mod) * avg;
			next_index = mod * (avg + 1) + (thread_id - mod + 1) * avg;
		}

		// Convert index into row and column numbers
		int l2x = index / 4;
		int l2y = index % 4;
		int l2xe = min((next_index - 1) / 4, 3);
		int l2ye = min((next_index - 1) % 4, 3);

		//----------------------------------------------------

		// box_list is a list each thread has, in the case of serial there is only 1 box_list
		// box_list contains the list of all the boxes, on each layer, that fmm must go through and calculate
		// This is done here so that fmm can just do calculations instead of figuring out in which
		// boxes to do the calculation
		vector < vector < pair < int, int > > > box_list;
		for(int i = 0; i <= depth_level; i++)
		{
			vector < pair<int, int > > dummy_v;
			box_list.push_back(dummy_v);
		}

		// This part wanders through the level 2 boxes this process got assigned to itself
		// It calls the function find_bottom_box(), this function goes to the bottom layer
		// children of this box and add them to this process's bottom_boxes list.
		int tcx = l2x;
		int tcy = l2y;
		for(int nn = next_index; index < nn; nn--)
		{
			if(tcy > 3)
			{
				tcy = 0;
				tcx++;
			}
			find_bottom_box(boxes[2][tcx][tcy], boxes, box_list, depth_level, tcx, tcy);
			tcy++;
		}

		#pragma omp barrier

		for( int step = 0; step < 1000; step++ )
		{
			#pragma omp barrier

			// Calls Fast Multipole Method
			omp_fmm(boxes, particles, 2, l2x, l2y, l2xe, l2ye);

			#pragma omp barrier

			//  move particles
			for(int i = 0; i < box_list[2].size(); i++)
			{
				int boxx = box_list[2][i].first;
				int boxy = box_list[2][i].second;
				for(int j = 0; j < boxes[2][boxx][boxy].element.size(); j++)
				{ 
					// Moves particle
					move(particles[boxes[2][boxx][boxy].element[j]]);	
					// Resets particle force
					particles[boxes[2][boxx][boxy].element[j]].ax = particles[boxes[2][boxx][boxy].element[j]].ay = 0;
				}
			}

			#pragma omp barrier

			// If thread is 0 send root particles to level 1
			if(thread_id == 0)
				sort_elements(world, particles, boxes, 0, 0);

			#pragma omp barrier

			// Send particles from level 1 to level 2
			if(thread_id < 4)
				sort_elements(boxes[1][x][y], particles, boxes, x, y);

			#pragma omp barrier

			// Send the particles to the rest of the boxes
			assign_elements(boxes, particles, 2, depth_level, l2x, l2y, l2xe, l2ye);



			//  save if necessary
			#pragma omp master
			if( fsave && (step%SAVEFREQ) == 0 )
				save( fsave, n, particles );
		}
	}

	simulation_time = read_timer( ) - simulation_time;

	printf("n = %d,\tsimulation time = %g seconds\n", n, simulation_time );

	free( particles );
	if( fsave )
		fclose( fsave );

	return 0;
}
